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Abstract 

Parallel MRI is a fast imaging technique that enables the acquisition of highly 
resolved images in space. It relies on fc-space undersampling and multiple re- 
ceiver coils with complementary sensitivity profiles in order to reconstruct a 
full Field-Of-View (FOV) image. The performance of parallel imaging mainly 
depends on the reconstruction algorithm, which can proceed either in the 
original fc-space (GRAPPA, SMASH) or in the image domain (SENSE-like 
methods). To improve the performance of the widely used SENSE algorithm, 
2D- or slice-specific regularization in the wavelet domain has been efficiently 
investigated. In this paper, we extend this approach using 3D- wavelet rep- 
resentations in order to handle all slices together and address reconstruction 
artifacts which propagate across adjacent slices. The extension also accounts 
for temporal correlations that exist between successive scans in functional 
MRI (fMRI). The proposed 4D reconstruction scheme is fully unsupervised 
in the sense that all regularization parameters are estimated in the maximum 
likelihood sense on a reference scan. The gain induced by such extensions 
is first illustrated on EPI image reconstruction but also measured in terms 
of statistical sensitivity during a fast event-related fMRI protocol. The pro- 
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posed 4D-UWR-SENSE algorithm outperforms the SENSE reconstruction at 
the subject and group-levels (15 subjects) for different contrasts of interest 
and using different parallel acceleration factors on 2 x 2 x 3mm^ EPI images. 

Keywords: pMRI, SENSE, frames, sparsity, fMRI, convex optimization, 
parallel computing, proximal methods 



1. Introduction 

Reducing scanning time in Magnetic Resonance Imaging (MRI) exams 
remains a worldwide challenging issue. The expected benefits of a faster 
acquisition can be summarized as follows: i) limit patient's exposure to the 
MRI environment either for safety or discomfort reasons, ii) maintain a strong 
robustness in the acquisition with respect to subject's motion artifacts, iii) 
limit geometric distortions or maintain high image quality and iv) acquire 
more spatially or temporally resolved images in the same or even reduced 



amount of time ||22L |36|]. The basic idea to make MRI acquisitions faster 
or to improve spatial resolution at fixed scanning time consists of reducing 
the amount of acquired samples in the fc-space and developing dedicated 
reconstruction pipelines. To achieve this goal, two main research avenues 
have been developed so far: i) parallel imaging that relies on a geometri- 
cal complementarity principle involving multiple receiver channel coils with 
complementary sensitivity profiles. This enables the reduction of the num- 
ber of fc-space lines to be acquired without degrading spatial resolution or 
truncating the Field of View (FOV), and thus requires the combination of 



reduced FOV coil-specific images to reconstruct the full FOV image [3J, [18 
ii) compressed sensing MRI that exploits the implicit sparsity in MR images 
to significantly undersample the /c-space and randomly select incoherent (or 
complementary) samples regarding their spectral contribution to the MR im- 
age [2^ . This approach which does not require any multiple channel coil and 
thus remains useable with birdcage ones will not be addressed in this work. 

The present paper is a contribution to parallel imaging at conventional 
magnetic field strength (e.g. 3 Tesla). At 3 Tesla, Parallel Magnetic Res- 
onance Imaging (pMRI) is mainly useful in reception, i.e. at the image 
reconstruction step using multiple channel coils, while higher magnetic fields 
(i.e. 7 Tesla) require parallel imaging at the transmission step as well |2ol. l3|. 

Many methods like GRAPPA (Generalized Autocalibrating Partially Par- 
allel Acquisitions) [l8| and SENSE (Sensitivity Encoding) [34] have been 
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proposed in the literature to reconstruct a full FOV image from multiple 
/c-space under-sampled images acquired on separate channels. The main dif- 
ference between these two classes of methods lies in the space on which they 
operate: GRAPPA performs multichannel full FOV reconstruction in the 
/c-space whereas SENSE carries out the unfolding process in the image do- 
main: all the under-sampled images are first reconstructed by inverse Fourier 
transform before combining them to unwrap the full FOV image. Another 
difference is that GRAPPA is autocalibrated, while SENSE needs a separate 
coil sensitivity estimation step based on a reference scan. Note however that 
an autocalibrated version of SENSE is also available for instance in Siemens 
scanners and called mSENSE hereafter. SENSE as well as GRAPPA methods 
may suflFer from strong artifacts when high values of acceleration factors R 
are considered in the imaging setup or when they are applied to Echo Planar 
Imaging (EPI), the sequence involved in fMRI experiments. These artifacts 
can drastically disturb subsequent statistical analysis such as brain activation 
detection. Regularized SENSE methods have been proposed in the literature 
to improve the robustness of the solution 23, 46, 4^,0,113]. Some of them ap- 
ply quadratic or Total Variation (TV) regularizations while others resort to 
regularization in the wavelet transform domain (e.g. UWR-SENSE [7]). The 
latter strategy has proved its eflSciency on the reconstruction of anatomical or 
functional (resting-state only) data jy, |7|. More recently, UWR-SENSE has 
been assessed on EPI images and compared with mSENSE on the same data 
acquired during a brain activation fMRI experiment [5]. This comparison 
was performed at the subject level on a few subjects only. 

Besides, most of the available reconstruction methods in the literature 
operate slice by slice and thus reconstruct each slice irrespective of its neigh- 
bors. Iterating over slices is thus necessary to get the whole 3D volume. This 
observation led us to consider 3D or whole brain reconstruction as a single 
step in which all slices are treated together by making use of 3D wavelet trans- 
forms and 3D sparsity promoting regularization terms in the wavelet domain. 
Following the same principle, an fMRI run usually consists of several hun- 
dred successive scans that are reconstructed independently one to another 
and thus implicitly assumed to be independent of each other. Iterating over 
all acquired 3D volumes remains the classical approach to reconstruct the 
4D or 3D + t dataset associated with an fMRI run. However, it has been 
shown for a long while that fMRI data are serially correlated in time even 
under the null hypothesis (i.e., ongoing activity only) (H, [13, 35]. To capture 
this dependence between successive time points, an autoregressive model has 
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demonstrated its relevance 4J, |45|, l3l| , especially when its parameters and 



optimal order can vary in space for instance with tissue type [43|, l30|]. Hence, 
it makes sense to account for this temporal structure at the reconstruction 
step. 

These two key ideas play a central role in the present paper by extend- 
ing the UWR-SENSE approach [7] through a fidelity data term combining 
all time points, which relies on a 3D wavelet transform and an additional 
regularization term along the temporal dimension of the 4D dataset in the 
image domain. The development of the proposed method (named 4D-UWR- 
SENSE) was made possible due to recent advances in nonsmooth convex 
optimization. Indeed, it is based on a Parallel ProXimal Algorithm (PPXA) 
which is different from the ones employed in 

4D-UWR-SENSE leads to improvements in the retrieval of a reliable 
Signal-to-Noise Ratio (SNR) between the acquired volumes and also in the 
enhancement of the detection of BOLD effects or evoked activations that oc- 
cur in response to the delivered stimuli in the fMRI experiment. The present 
paper therefore aims at demonstrating that the 4D-UWR-SENSE approach 
outperfoms its SENSE-like alternatives not only in terms of reduced artifacts, 
but also in terms of statistical sensitivity at the voxel and cluster levels, in 
intra-subject and group studies. The rest of this paper is organized as fol- 
lows. Section [2] recalls the general parallel MRI framework. We then describe 
the proposed reconstruction algorithm in Section [3] before illustrating related 
experimental results in Section [H Finally, some discussions and conclusions 
are drawn in Section [51 
2. Parallel imaging in MRI 

In parallel MRI, an array of L coils is employed tojneasure the spin 
density p into the object under investigatioij^. The signal di received by each 
coil ^ (1 < ^ < L) is the Fourier transform of the desired 2D field p G M^^^ 
on the specified FOV weighted by the coil sensitivity profile 5^, evaluated at 
some location kr = {kx^ kyY in the /c-space: 

UK) = I p(r)s,(r)e-*2xfe> ^ ~ ^^^^^ 



where n£{kr) is a coil-dependent additive zero-mean Gaussian noise, which is 
independent and identically distributed (iid) in the /c-space, and r = (x, yY G 



The overbar is used to distinguish the "true" data from a generic variable. 
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X X y is the spatial position in the image domain {-^ being the transpose 
operator). The size of the reduced FOV acquired data in the fc-space 
clearly depends on the sampling scheme. In contrast with cardiac imaging 
where the heart motion is significant during scanning, a Cartesian coordinate 
system is generally adopted in the neuroimaging context. In parallel MRI, 
the sampling period along the phase encoding direction is R times larger 
than the one used for conventional acquisition, R < L being the reduction 
factor. To recover full FOV images, many algorithms have been proposed 
but only SENSE-hke H and GRAPPA-hke |18| methods are provided by 
scanner manufacturers. In what follows, we focus on SENSE-like methods 
operating in the spatial domain. 

Let ^ ^ be the aliasing period and y the position in the image 
domain along the phase encoding direction. Let x be the position in the 
image domain along the frequency encoding direction. A 2D inverse Fourier 
transform allows us to recover the measured signal in the spatial domain. 
By accounting for the fc-space undersampling at i?-rate, the inverse Fourier 
transform gives us the spatial counterpart of Eq. ([T]) in matrix form: 

d{r) = S{r)p{r) + n{r), (2) 

where 

ni{x,y) 



Sir) 



si{x,y) ... si{x,y + (R- l)Ay) 



SL{x,y) ... SL{x,y + [R- l)Ay) 



n[r] 



-( \ ^ 
p{r) = 



p{x,y) 

p{x, y + Ay) 
p{x,y + {R-l)Ay) 



and d{r) 



di{x,y) 
d2ix,y) 

dLix,y) 



(3) 



Based upon this model, the reconstruction step consists of solving Eq. ([2j) 
so as to recover p(r) from d{r) and an estimate of S{r) at each spatial posi- 
tion r = (x, yy. The spatial mixture or sensitivity matrix S{r) is estimated 



using a reference scan and varies according to the coil geometry. Note that 
the coil images (rf^)i<KL as well as the sought image p are complex- valued, 
although \p\ is only considered for visualization purposes. The next section 
describes the widely used SENSE algorithm as well as its regularized exten- 
sions. 

3. Reconstruction algorithms 

3.1. ID-SENSE 

In its simplest form, SENSE imaging amounts to solving a one-dimensional 
inversion problem due to the separability of the Fourier transform. Note 
however that this inverse problem admits a two-dimensional extension in 
3D imaging sequences like Echo Volume Imaging (EVI) |36|] where under- 
sampling occurs in two fc-space directions. The ID-SENSE reconstruction 
method [33] actually minimizes a Weighted Least Squares (WLS) criterion 
JwLS given by: 

JwLsip) = E II - S(^)P(^) III-' (4) 

rG{l,...,X}x{l,...,y/i?} 



where || • ||^-i = \/ (•)^^~''"(-), and the noise covariance matrix ^ is usually 
estimated based on L acquired images (d^)i<^<L from all coils without radio 
frequency pulse. Hence, the SENSE full FOV image is nothing but the max- 
imum likelihood estimate under Gaussian noise assumption, which admits 
the following closed- form expression at each spatial position r: 

PwLS W = {S^{r)^-'S{r)y (r)^-' d{r) , (5) 

where (•)^ (resp. (•)^) stands for the transposed complex conjugate (resp. 
pseudo-inverse). It should be noticed here that the described ID-SENSE re- 
construction method has been designed to reconstruct one slice (2D image). 
To reconstruct a full volume, the ID-SENSE reconstruction algorithm has to 
be iterated over all slices. 

In practice, the performance of the ID-SENSE method is limited because of i) 
the presence of distortions in the measurements d(r), ii) the ill-conditioning 
of S'(r), in particular at locations r close to the image center and in) the 
presence of errors in the estimation of S{r) mainly at brain/air interfaces. To 
enhance the robustness of the solution to this ill-posed problem, a regulariza- 
tion is usually introduced in the reconstruction process. To improve results 
obtained with quadratic regularization techniques [2^, 46], edge-preserving 



6 



regularization has been widely investigated in the pMRI reconstruction hter- 
ature. For instance, reconstruction methods based on Total Variation (TV) 



regularization have been proposed in a number of recent works like [2l|, |25| . 
However, TV is mostly adapted to piecewise constant images, which are not 
always acurate models in MRI, especially in fMRI. As investigated by Chaari 
et al I?! and Liu et al [24], regularization in the Wavelet Transform (WT) 
domain is a powerful tool to improve SENSE reconstruction. In what follows, 
we summarize the principles of the wavelet-based regularization approach. 

3.2. Proposed wavelet Regularized- SENSE 

Akin to [7] where a regularized reconstruction algorithm relying on 2D 
separable WTs was investigated, to the best of our knowledge, all the existing 
approaches in the pMRI regularization literature proceed slice by slice. The 
drawback of this strategy is that no spatial continuity between adjacent slices 
is taken into account since the slices are processed independently. Moreover, 
since the whole brain volume has to be acquired several times in an fMRI 
study, iterating over all the acquired 3D volumes is then necessary in order 
to reconstruct a 4D data volume corresponding to an fMRI session. Conse- 
quently, the 3D volumes are supposed independent whereas fMRI time-series 
are serially correlated in time because of two distinct effects: the BOLD sig- 
nal itself is a low-pass filtered version of the neural activity, and physiological 
artifacts make the fMRI time points strongly dependent. For these reasons, 
modeling temporal dependence across scans at the reconstruction step may 
impact subsequent statistical analysis. This has motivated the extension of 
the wavelet regularized reconstruction approach in 0] in order to: 

• account for 3D spatial dependencies between adjacent slices by using 
3D WTs, 

• exploit the temporal dependency between acquired 3D volumes by ap- 
plying an additional regularization term along the temporal dimension 
of the 4D dataset. 

This additional regularization will help in increasing the Signal to Noise Ra- 
tio (SNR) through the acquired volumes, and therefore enhance the reliability 
of the statistical analysis in fMRI. These temporal dependencies have also 
been used in the dynamic MRI literature in order to improve the reconstruc- 
tion quality in conventional MRI [39]. However, since the imaged object 
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geometry in the latter context generally changes during the acquisition, tak- 
ing into account the temporal regularization in the reconstruction process is 
very difficult. 

To deal with a 4D reconstruction of the A^^ acquired volumes, we will first 
rewrite the observation model in Eq. ([2j) as follows: 

d\r) = S{r)p\r) + n\r), (6) 

where t G {1,...,A^^} is the acquisition time and r = {x^y^z) is the 3D 
spatial position, z G {1, . . . , Z} being the position along the third direction 
(slice selection one). 

At a given time t, the full FOV 3D complex- valued image of size XxY xZ 
can be seen as an element of the Euclidean space with K = X xY x Z 
endowed with the standard inner product (•!•) and norm || • ||. We em- 
ploy a dyadic 3D orthonormal wavelet decomposition operator T over jmax 
resolution levels. The coefficient field resulting from the wavelet decompo- 
sition of a target image is defined as = (CL (Co j)oGO,i<j<jmax) with 
o e O = {0,1}3\ {(0,0,0)}, C* = (CU)i»,_ and ="(C|Ji<it<K, 
where Kj = K2 is the number of wavelet coefficients in a given subband at 
resolution j (by assuming that X, Y and Z are multiple of 2-^"'^^). Adopting 
such a notation, the wavelet coefficients have been reindexed so that de- 
notes the approximation coefficient vector at the resolution level jmax, while 
denotes the detail coefficient vector at the orientation o and resolution 
level j. Using 3D dyadic WTs allows us to smooth reconstruction artifacts 
along the slice selection direction, which is not possible using a slice by slice 
operating approach. 

The proposed regularization procedure relies on the introduction of two 
penalty terms. The first penalty term describes the prior spatial knowl- 
edge about the wavelet coefficients of the target solution and it is expressed 
as: 

Nr -^Jmax Jmax 

9(0 = E [ E + E E E '^oAC,j,k)] , (7) 

t=i k=i oeo j=i k=i 

where ( = ((^^, (^^, . . . , C^^) and we have, for every o G O and j G {1, . . . , jmax}, 
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where $^^(0 = <^|Re(^-/i,,,-)l + TlRe(^-//o,i)P and $^-(0 = <|Im(e- 

Aio,,)| + T|Im(e - with /i,,, = + G C, and /5,^|, a^-, 

are some positive real constants. Hereabove, Re(-) and Im(-) (or -^^ and 
•^"^) stand for the real and imaginary parts, respectively. A similar model is 
adopted for the approximation coefficients. The second regularization term 
penalizes the temporal variation between successive 3D volumes: 

Nr 

t=2 

where T* is the 3D wavelet reconstruction operator. The prior parameters 

= «,«^^), /3o,j = {P^^j.P'^). f^o,j = (/^^^^o^S)' « e M+ and p e 
[l,+oc[ are unknown and they need to be estimated 
(see [Appendix A\ ) . 

The operator T* is then applied to each component of ( to obtain the 
reconstructed 3D volume related to the acquisition time t. It should be 
noticed here that other choices for the penalty functions are also possible pro- 
vided that the convexity of the resulting optimality criterion is ensured. This 
condition enables the use of fast and efficient convex optimization algorithms. 
Adopting this formulation, the minimization process plays a prominent role 
in the reconstruction process, as detailed in [Appendix B[ 

4. Experimental validation in fMRI 

This section is dedicated to the experimental validation of the 4D-UWR- 
SENSE reconstruction algorithm we proposed in Section 13. 2[ Results of 
subject and group-level fMRI statistical analyses are compared for two re- 
construction pipelines: one available on the Siemens workstation and our 
own pipeline involving for the sake of completeness either the early UWR- 
SENSE [7] or the 4D-UWR-SENSE version of the proposed pMRI recon- 
struction algorithm. In what follows, we first describe the fMRI acquisition 
setup and the experimental design. 

4.1. Experimental data 

For validation purpose, we acquired fMRI data on a 3 T Siemens Trio 
magnet using a Gradient-Echo EPI (GE-EPI) sequence {TE = 30 ms, TR = 
2400 ms, slice thickness = 3 mm, transversal orientation, FOV = 192 mm^) 
during a cognitive localizer |[32|] protocol. This experiment has been designed 
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to map auditory, visual and motor brain functions as well as higher cognitive 
tasks such as number processing and language comprehension (listening and 
reading). It consisted of a single session of A^^ = 128 scans. The paradigm was 
a fast event-related design comprising sixty auditory, visual and motor stim- 
uli, defined in ten experimental conditions (auditory and visual sentences, 
auditory and visual calculations, left/right auditory and visual clicks, hori- 
zontal and vertical checkerboards). An L = 32 channel coil was used to en- 
able parallel imaging. Ethics approval was given by the local research ethics 
committee, and fifteen subjects gave written informed consent for participa- 
tion. For each subject, fMRI data were collected at the 2x2 mm^ spatial 
in-plane resolution using different reduction factors (i? = 2 or i? = 4). Based 
on the raw data files delivered by the scanner, reduced FOV EPI images were 
reconstructed as detailed in Fig. [H This reconstruction requires two specific 
treatments: 

i) k- space regridding to account for the non-uniform /c-space sampling 
during readout gradient ramp, which occurs in fast MRI sequences like 
GE-EPI; 

ii) Nyquist ghosting correction to remove the odd-even echo inconsistencies 
during fc-space acquisition of EPI images. 



Siemens MRI 


Reading raw FID 


scanner 


data 





Regridding 




Nyquist ghosting 




IFFT 








correction 







Reduced FOV 



images 

Figure 1: Reconstruction pipeline of reduced FOV EPI images from the raw FID data. 



Once the reduced FOV images are available, the proposed pMRI 4D-UWR- 
SENSE algorithm and its early UWR-SENSE version have been utilized in 
a final step to reconstruct the full FOV EPI images and compared to the 
hiSENSeI Siemens solution. For the wavelet-based regularization, dyadic (M = 
2) Symmlet orthonormal wavelet bases fl5| associated with filters of length 
8 have been used over jmax = 3 resolution levels. The reconstructed EPI 
images then enter in our fMRI study in order to measure the impact of the 
reconstructor choice on brain activity detection. Note also that the pro- 
posed reconstruction algorithm requires the estimation of the coil sensitivity 



^ SENSE reconstruction implemented by the Siemens scanner 
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maps (matrix S{-) in Eq. ([2j)). As proposed in j3J], the latter were esti- 
mated by dividing the coil-specific images by the module of the Sum Of 
Squares (SOS) images, which are computed from the specific acquisition of 
the /c-space center (24 lines) before the A^^ scans. 

Fig. [2] compares the two pMRI reconstruction algorithms to illustrate 
on axial, coronal and sagittal slices how the mSENSE reconstruction artifacts 
have been removed using the 4D-UWR-SENSE approach. The mSENSE re- 
constructed images acutally present large artifacts located both at the center 
and boundaries of the brain in sensory and cognitive regions (temporal lobes, 
frontal and motor cortices, ...). This results in SNR loss and thus may have 
a dramatic impact for activation detection in these brain regions. Note that 
these conclusions are reproducible across subjects although the artifacts may 
appear on different slices (see red circles in Figjjj). It is also worth men- 
tioning that in contrast to the Siemens reconstructor our pipeline does not 
involve any spatial filtering step to improve signal homogeneity across the 
brain: this explains why the images shown in Fig. [2] and delivered by our 
4D-UWR-SENSE algorithm seem less homogeneous. However, bias field cor- 
rection can be applied if necessary using specific tools such as those available 
in BrainVISAl. 

Regarding the computation burden, the mSENSE algorithm is carried out 
on-line and remains compatible with real time processing. On the other hand, 
our pipeline is carried out off-line and requires much more computation. For 
illustration purpose, on a biprocessor quadicore Intel Xeon CPU® 2.67GHz, 
one EPI slice is reconstructed in 4 s using the UWR-SENSE algorithm. Using 
parallel computing strategy and multithreading (through the OMP library), 
each EPI volume consisting of 40 slices is reconstruced in 22 s. This makes 
the whole series of 128 EPI images available in about 47 min. In contrast, 
the proposed 4D-UWR-SENSE achieves the reconstruction of the series in 
about 40 min, but requires larger memory space due to large data volume 
processed simultaneously. 

4.2. fMRI data pre-processings 

Irrespective of the reconstruction pipeline, the full FOV fMRI images 
were then preprocessed using the SPM5 softwar^: preprocessing involves 



"^http : / /brainvisa . info 

^ http : // ww w . f il . ion . ucl . ac . uk/spm/sof tware/spin5/| 
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R = 2 



R = 4 



mSENSE 4D-UWR-SENSE 



Axial 



Coronal 



Sagittal 



Axial 



Coronal 






Sagittal 



Figure 2: Axial, Coronal and Sagittal reconstructed slices using mSENSE and 4D-UWR- 
SENSE for R = 2 and R = 4 with 2x2 mm^ in-plane spatial resolution. Red circles and 
ellipsoids indicate the position of reconstruction artifacts using mSENSE. 



realignment, correction for motion and differences in slice acquisition time, 
spatial normalization, and smoothing with an isotropic Gaussian kernel of 
4mm full-width at half-maximum. Anatomical normalization to MNI space 
was performed by coregistration of the functional images with the anatomical 
Tl scan acquired with the thirty two channel-head coil. Parameters for the 
normalization to MNI space were estimated by normalizing this scan to the 
Tl MNI template provided by SPM5, and were subsequently applied to all 
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functional images. 



4-3. Subject-level analysis 

A General Linear Model (GLM) was constructed to capture stimulus- 
related BOLD response. As shown in Fig. [3l the design matrix relies on 
ten experimental conditions and thus made up of twenty one regressors cor- 
responding to stick functions convolved with the canonical Haemodynamic 
Response Function (HRF) and its first temporal derivative, the last regressor 
modelling the baseline. This GLM was then fitted to the same acquired im- 
ages but reconstructed using either the Siemens reconstructor or our own 
pipeline, which in the following is derived from the early UWR-SENSE 
method [7] and from its 4D-UWR-SENSE extension we propose here. 




10 15 



Figure 3: (a): design matrix and the Lc-Rc contrast involving two conditions (grouping 
auditory and visual modalities); (b): design matrix and the aC-aS contrast involving four 
conditions (sentence, computation, left click, right click). 



Here, contrast estimate images for motor responses and higher cognitive 
functions (computation, language) were subjected to further analyses at the 
subject and group levels. These two contrast are complementary since the ex- 
pected activations lie in diflFerent brain regions and thus can be diflFerentially 
corrupted by reconstruction artifacts as outlined in Fig. [2l More precisely, 
we studied: 



• the Auditory computation vs. Auditory sentence (aC-aS) con- 
trast which is supposed to elicit evoked activity in the frontal and pari- 
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etal lobes, since solving mental arithmetic task involves working mem- 
ory and more specifically the intra-parietal sulcus [16]: see Fig. [3](b); 

• the Left click vs. Right click (Lc-Rc) contrast for which we expect 
evoked activity in the right motor cortex (precentral gyrus, middle 
frontal gyrus). Indeed, the Lc-Rc contrast defines a compound com- 
parison which involves two motor stimuli which are presented either 
in the visual or auditory modality. This comparison aims therefore at 
detecting lateralization effect in the motor cortex: see Fig. [3](a). 

Interestingly, these two contrasts were chosen because they summarized 
well different situations (large vs small activation clusters, distributed vs fo- 
cal activation pattern, bilateral vs unilateral activity) that occurred for this 
paradigm when looking at sensory areas (visual, auditory, motor) or regions 
involved in higher cognitive functions (reading, calculation). In the follow- 
ing, our results are reported in terms of Student-t maps thresholded at a 
cluster-level p = 0.05 corrected for multiple comparisons according to the 
Family Wise Error Rate (FWER) [2i, 1^. Complementary statistical tables 
provide corrected cluster and voxel-level values, maximal t-scores and cor- 
responding peak positions both for i? = 2 and i? = 4. Note that clusters are 
listed in a decreasing order of significance. 

Concerning the aC-aS contrast. Fig. Hftop] shows for the most signifi- 
cant slice and R = 2 that all pMRI reconstruction algorithms succeed in 
finding evoked activity in the left parietal and frontal cortices, more pre- 
cisely in the inferior parietal lobule and middle frontal gyrus according to 
the AAL templat^. Table [1] also confirms a bilateral activity pattern in 
parietal regions for R = 2. Moreover, for i? = 4 Fig. [Hbottom] illustrates 
that our pipehne (UWR-SENSE and 4D-UWR-SENSE) and preferentially 
the proposed 4D-UWR-SENSE scheme enables to retrieve reliable frontal 
activity elicited by mental calculation, which is lost by the the mSENSE al- 
gorithm. From a quantitative viewpoint, the proposed 4D-UWR-SENSE 
algorithm finds larger clusters whose local maxima are more significant than 
the ones obtained using mSENSE and UWR-SENSE, as reported in Table [H 
Concerning the most significant cluster for i? = 2, the peak positions re- 
main stable whatever the reconstruction algorithm. However, examining 
their significance level, one can first measure the benefits of wavelet-based 



^available in the xjView toolbox of SPM5. 
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regularization when comparing UWR-SENSE with mSENSE results and then 
additional positive effects of temporal regularization and 3D wavelet decom- 
position when looking at the 4D-UWR-SENSE results. These benefits are 
also demonstrated for i? = 4. 
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Figure 4: Subject-level student-t maps superimposed to anatomical MRI for the aC-aS 
contrast. Data have been reconstructed using the mSENSE, UWR-SENSE and 4D-UWR- 
SENSE, respectively. Neurological convention: left is left. The blue cross shows the 
maximum activation peak. 



Fig. [5] illustrates another property of the proposed pMRI pipeline, ie 
its robustness to the between-subject variability. Indeed, when comparing 
subject-level student-t maps reconstructed using the different pipelines {R = 
2), it can be observed that the mSENSE algorithm fails to detect any activation 
cluster in the expected regions for the second subject (see Fig. [5[bottom]). 
In contrast, our 4D-UWR-SENSE method retrieves more coherent activity 
while not exactly at the same position as for the first subject. 

For the Lc-Rc contrast on the data acquired with i? = 2, Fig. [6[top] 
shows that all reconstruction methods enable to retrieve expected activation 
in the right precentral gyrus. However, when looking more carefully at the 
statistical results (see Table [2j) , our pipeline and more preferentially the 4D- 
UWR-SENSE algorithm retrieves an additional cluster in the right middle 
frontal gyrus. On data acquired with i? = 4, the same Lc-Rc contrast elicits 
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Table 1: Significant statistical results at the subject-level for the aC-aS contrast (corrected 
for multiple comparisons at p = 0.05). Images were reconstructed using the mSENSE, UWR- 
SENSE and 4D-UWR-SENSE algorithm for R = 2 and R = 4. 
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similar activations, ie in the same region. As demonstrated in Fig.[6[bottom], 
this activity is enhanced when pMRI reconstruction is performed with our 
pipeline. Quantitative results in Table [2] confirms numerically what can be 
observed in Fig.O larger clusters with higher local t-scores are detected using 
the 4D-UWR-SENSE algorithm, both for i? = 2 and R = 4. Also, a larger 
number of clusters is retrieved for i? = 2 using wavelet-based regularization. 

Fig. [7] reports on the robustness of the proposed pMRI pipeline to the 
between-subject variability for this motor contrast. Since sensory functions 
are expected to generate larger BOLD effects (higher SNR) and appears more 
stable, our comparison takes place at i? = 4. Two subject-level student-t 
maps reconstructed using the different pMRI algorithms are compared: in 
Fig. [71 For the second subject, one can observe that the mSENSE algorithm 
fails to detect any activation cluster in the right motor cortex. In contrast. 
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Figure 5: Between-subject variability of detected activation for the aC-aS contrast at 
R = 2. Neurological convention. The blue cross shows the activation peak. 



Table 2: Significant statistical results at the subject-level for the Lc-Rc contrast (corrected 
for multiple comparisons at p = 0.05). Images were reconstructed using the mSENSE, UWR- 
SENSE and 4D-UWR-SENSE algorithms for R = 2 and = 4. 
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our 4D-UWR-SENSE method retrieves more coherent activity for this second 
subject in the expected region. 

To summarize, on these two contrasts our 4D-UWR-SENSE algorithm 
always outperforms the alternative reconstruction methods in terms of sta- 
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Figure 6: Subject-level student-t maps superimposed to anatomical MRI for the Lc-Rc 
contrast. Data have been reconstructed using the mSENSE, UWR-SENSE and 4D-UWR- 
SENSE, respectively. Neurological convention. The blue cross shows the activation peak. 



tistical significance (number of clusters, cluster extent, peak values,...) but 
also in terms of robustness. 

4.4. Group-level analysis 

Due to between-subject anatomical and functional variability, group-level 
analysis is necessary in order to derive robust and reproducible conclusions 
at the population level. For this validation, random effect analyses (RFX) 
involving fifteen healthy subjects have been conducted on the contrast maps 
we previously investigated at the subject level. More precisely, one-sample 
Student-t test was performed on the subject-level contrast images (eg, Lc-Rc, 
aC-aS,... images) using SPM5. 

For the aC-aS contrast. Maximum Intensity Projection (MIP) student- 
t maps are shown in Fig. [8l First, they illustrate that irrespective of the 
reconstruction method larger and more significant activations are found on 
datasets acquired with R = 2 given the better SNR. Second, for i? = 2, 
visual inspection of Fig. [8[top] confirms that only the 4D-UWR-SENSE al- 
gorithm allows us to retrieve significant bilateral activations in the parietal 
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Figure 7: Between-subject variability of detected activation for the Lc-Rc contrast at 
R = 4. Neurological convention. The blue cross shows the activation peak. 



cortices (see axial MIP slices) in addition to larger cluster extent and a gain 
in significance level for the stable clusters across the different reconstructors. 
Similar conclusions can be drawn when looking at Fig. [8[bottom] for i? = 4. 
Complementary results are available in Table [3] for i? = 2 and i? = 4 and 
numerically confirms this visual comparison: 

• Whatever the reconstruction method in use, the statistical performance 
is much more significant using i? = 2, especially at the cluster level since 
the cluster extent decreases by one order of magnitude. 

• Voxel and cluster-level results are enhanced using the 4D-UWR-SENSE 
approach instead of the mSENSE reconstruction or its early UWR-SENSE 
version. 

Fig. M reports similar group-level MIP results for i? = 2 and i? = 4 con- 
cerning the Lc-Rc contrast. It is shown that whatever the acceleration factor 
R in use, our pipeline enables to detect a much more spatially extended ac- 
tivation area in the motor cortex. This visual inspection is quantitatively 
confirmed in Table [4] when comparing the detected clusters using our 4D- 
UWR-SENSE approach with those found by mSENSE, again irrespective of R. 
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Figure 8: Group-level student-t maps for the aC-aS contrast where data have been recon- 
structed using the mSENSE, UWR-SENSE and 4D-UWR-SENSE for = 2 and R = 4. 
Neurological convention. Red arrows indicate the global maximum activation peak. 



Finally, the 4D-UWR-SENSE algorithm outperforms the UWR-SENSE one, 
which corroborates the benefits of the proposed spatio-temporal regulariza- 
tion scheme. 

5. Discussion and conclusion 

The contribution of the present paper was twofold. First, we proposed 
a novel reconstruction method that relies on 3D wavelet transform and ac- 
counts for temporal dependencies in successive fMRI volumes. Second, our 
particular interest was to demonstrate that when artifacts are superimposed 
to brain activation, this directly impacts subsequent brain activity detection. 
In this context, we showed that the choice of the parallel imaging reconstruc- 
tion algorithm impacts the statistical sensitivity in fMRI at the subject and 
group-levels and may enable whole brain neuroscience studies at high spatial 
resolution. 

Practically speaking, we showed that whole brain acquisition can be rou- 
tinely used at a spatial in-plane resolution of 2 x 2mm^ in a short and constant 
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Table 3: Significant statistical results at the group-level for the aC-aS contrast (corrected 
for multiple comparisons at p = 0.05). Images were reconstructed using the mSENSE, 
UWR-SENSE and 4D-UWR-SENSE algorithms for R = 2 and R = 4. 
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repetition time (TR = 2.4s) provided that a reliable pMRI reconstruction 
pipeline is chosen. In this paper, we demonstrated that our 4D-UWR-SENSE 
reconstruction algorithm meets this desired property. To draw this conclu- 
sion, our comparison took place at the statistical analysis level and relied 
on quantitative criteria (voxel- and cluster-level corrected p- values, t-scores, 
peak positions) both at the subject and group levels. In particular, we showed 
that our 4D-UWR-SENSE contribution outperforms both its UWR-SENSE 
ancestor Q and the Siemens mSENSE reconstruction in terms of statistical 
significance and robustness. This emphasized the benefits of combining tem- 
poral and 3D regularizations in the wavelet domain. Interestingly, we exhib- 
ited a most significant gain in the more degraded situation {R = 4) due to 
the positive impact of regularization. The strength of our conclusions lies 
in the reasonable size of our datasets since the same cohort participated to 
parallel imaging acquisitions using two different acceleration factors {R = 2 
and R 4). 

At this spatio-temporal compromise (2 x 2 x 3mm^ and TR = 2.4 s). 
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Figure 9: Group-level student-t maps for the Lc-Rc contrast where data have been recon- 
structed using the mSENSE, UWR-SENSE and 4D-UWR-SENSE for = 2 and R = A. 
Neurological convention. Red arrows indicate the global maximum activation peak. 



Table 4: Significant statistical results at the group-level for the Lc-Rc contrast (corrected 
for multiple comparisons at p = 0.05). Images were reconstructed using the mSENSE, 
UWR-SENSE and 4D-UWR-SENSE algorithms for = 2 and = 4. 
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we also illustrated the impact of increasing the acceleration factor (passing 
from i? = 2 to i? = 4) on the statistical sensitivity at the subject and group 
levels for a given reconstruction algorithm. We performed this comparison 
to anticipate what could be the statistical performance for detecting evoked 
brain activity on data requiring this acceleration level, such as high spatial 
resolution EPI images (eg 1.5x 1.5mm^ as in-plane resolution) acquired at the 
same short TR. Our conclusions were balanced depending on the contrast of 
interest: when looking at the aC-aS contrast, which involves in the fronto- 
parietal circuit, it turned out that i? = 4 was not reliable enough to recover 
group-level significant activity at 3 Tesla: the SNR loss was too important 
and should be compensated by an increase of the static magnetic field (eg 
passing from 3 to 7 Tesla). However, the situation is less dramatic or even 
acceptable for the Lc-Rc motor contrast, which elicits activation in motor 
regions: Our results brought evidence that the 4D-UWR-SENSE approach 
enabled the use of i? = 4 for this contrast. 

To summarize, the compromise between the acceleration factor and spa- 
tial in-plane resolution should be selected with care depending on the regions 
involved in the fMRI paradigm. As a consequence, high resolution fMRI 
studies can be conducted using high speed acquisition (short TR and large 
R value) provided that the expected BOLD eflFect is strong, as experienced 
in primary motor, visual and auditory cortices. Of course, the use of an op- 
timized reconstruction method such as the one proposed is a pre-requisite to 
shift this compromise towards larger R values and higher spatial resolution 
and could be optimally combined with ultra high magnetic fields. 

A direct extension of the present work, which is actually in progress, 
consists in studying the impact of tight frames instead of wavelet basis to de- 
fine shift-invariant transformation |7|]. However, unsupervised reconstruction 
becomes more challenging in this framework since the estimation of hyper- 
parameters becomes cumbersome (see [8] for details). Ongoing work will 
concern the combination of the present contribution with the joint detec- 
tion estimation approach of evoked activity 0, 0] to go beyond the GLM 
framework and measure how the pMRI reconstruction algorithm also impacts 
HRF estimation. Another extension would concern the combination of our 
wavelet-regularized reconstruction with the WSPM approach [4ol] in which 
statistical analysis is directly performed in the wavelet transform domain. 
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Appendix A. Maximum likelihood estimation of regularization pa- 
rameters 



A rigorous way of addressing the regularization parameter choice would 
be to consider that the sum of the regularization functions g and h cor- 
responds to the minus-log-likelihood of a prior distribution /(•;©) where 

imize the integrated likelihood of the data. This would however entail two 
main difficulties. On the one hand, this would require to integrate out the 
sought image decomposition ( and to iterate between image reconstruction 
and hyper-parameter estimation. Methods allowing us to perform this task 
are computationally intensive [17]. On the second hand, the partition func- 
tion of the distribution /(•;©) does not take a closed form and we would 



thus need to resort to numerical methods |4lL 1381. l37l| to compute it. To alle- 



viate the computational burden, akin to ll9|] we shall proceed differently by 
assuming that a^reference full FOV image p is available, and so is its wavelet 
decomposition ( = Tp. In practice, our reference image p is obtained using 
ID-SENSE reconstruction at the same R value. We then apply an approxi- 
mate ML procedure which consists of estimating separately the spatial and 
temporal parameters. Although this approach is not optimal from a theoret- 
ical standpoint, it is quite simple and it was observed to provide satisf actor 
results in practice. Alternative solutions based on Monte Carlo methods 



or Stein's principle llOj] can also be thought of, at the expense of an additional 



computational complexity. 

Appendix A.l. Spatial regularization parameters 

For the spatial hyper-parameter estimation task, we will assume that the 
real and imaginary parts of the wavelet coefficient^ are modelled by the 
following Generalized Gauss-Laplace (GGL) distribution: 




V^^M, /(e;/x,a,/3) = ^/^ -^-^^ . (A.l) 



For each resolution level j and orientation o, /i^^, afj and /3^j are estimated 

from (^^ j as follows (we proceed similarly to estimate and Pl"^ by 

replacing Re(-) by Im(-)): 



similar approach is adopted for the approximation coefficients. 
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= argmax /(Re(C,,); ^, a, /?) 

= argmax V log /(Re(Co,j,fc); yU, «, /?) 

(M,a,/3)eKxR+xR!j_ 

= arg min ^ « V \^e{Co,j,k - lA |+x |Re(Co,i,fc - ^) P 

This three-dimensional minimization problem does not admit a closed form 
solution. Hence, we can compute the ML estimated parameters using the 
zero-order Powell optimization method [l2l]. 

Appendix A. 2. Temporal regularization parameter 

For the temporal hyper-parameter estimation task, we will assume that, 
at a given voxel, the temporal noise is distributed according to the following 
generalized Gaussian (GG) distribution: 

l/p -K\e\P 

Akin to the spatial hyper-parameter estimation, reference images {p^)i<t<Nr 
are made available based on a ID-SENSE reconstruction, where 
Vt G {1, . . . , A^^}, = T*('^ We consider that at spatial position r, the tem- 
poral noise vector = [p^(r) — p^(r), p^(r) — p^(r), . . . , p-^^(r) — p^'^~^{r)Y 
is a realization of a full independent GG prior distribution and we adjust 
the temporal hyper-parameter vector directly from it. It should be 

noted here that the considered model for the temporal noise accounts for 
correlations between successive observations usually considered in the fMRI 
literature. It also presents more flexibility than the Gaussian model, which 
corresponds to the particular case when p = 2. Estimates 'k and p of the 
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parameters are then obtained as follows: 



(k3)= argmax f{er]i^,p) 

(/t,p)GM+x[l,+cx)[ 

= argmax log /(e^; 

(^,p)gM+x[1,+cx)[ 

= argmin J] |p*+i(r) - p*(r)r - (iV, - 1) log (-^^) . 

Kp)6R+x[i,+oo[ \2T{l/p)J 

(A.4) 

Note that in the above minimization, for a given value of the optimal value 
of K admits the following closed form: 

A zero-order Powell optimization method can then be used to solve the re- 
sulting one-variable minimization problem. To reduce the computational 
complexity of this estimation, it is only performed on the brain mask, and 
the temporal regularization parameter k is set to zero for voxels belonging 
to the image background. 

Appendix B. Optimization procedure for the 4D reconstruction 

We first recall that 

:^st(C) = Jtwls(C) + ^(C) + h{C) (B.l) 
where JTtwls is defined as 

Nr 

t=l 

Nr 

= E E \\d^{r)-S{r){T*C){r)\\l-.. 

t=l rG{l,...,X}x{l,...,y/i?}x{l,...,Z} 

(B.2) 

The minimization of ^st is performed by resorting to the concept of proxim- 



ity operators [28[] , which was found to be fruitful in a number of recent works 
in convex optimization 0, Q, 13]. In what follows, we recall the definition 
of a proximity operator: 
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Definition Appendix B.l [28|] Let To{x) be the class of lower semicon- 
tinuous convex functions from a separable real Hilbert space x to] — oo^ +oc] 
and let ip G ro(x)- For every x G X; the function + \\ • — x|p/2 achieves its 
infimum at a unique point denoted by prox^x. The operator prox^ ' X ^ X 
is the proximity operator of Lp. 

In this work, as the observed data are complex- valued, the definition of 
proximity operators is extended to a class of convex functions defined for 
complex- valued variables. For the function 

^] -oc,+oc] (B.3) 
X ^ (t)^^{Re{x)) + (/)^^(Im(x)), 

where and 0^"^ are functions in ro(M^) and Re(x) (resp. Im(x)) is the 
vector of the real parts (resp. imaginary parts) of the components of x G C^, 
the proximity operator is defined as 

prox^: C^^C^ (B.4) 
X ^ prox^Re(Re(x)) + zprox^im(Im(x)). 

Let us now provide the expression of the proximity operators involved in our 
reconstruction problem. 

Appendix B.l. Proximity operator of the data fidelity term 

According to standard rules on the calculation of proximity operators flil . 
Table 1.1], the proximity operator of the data fidelity term Jwls is given for 
every vector of coefficients (with t G {1, . . . , A^}) by proxj^^g(C^) = Tu^^ 
where the image is such that Vr G {1, . . . , X} x {1, . . . , y/^} x {1, . . . , Z}, 

^^(r) = (/^ + 25^(r)*-i5(r))"'(p^(r) + 25^(r)*-id^(r)), (B.5) 

where = T*^^ 

Appendix B.2. Proximity operator of the spatial regularization function 

According to [7], for every resolution level j and orientation o, the prox- 
imity operator of the spatial regularization function is given by 



+ ^ ''^""^y^ max{|Im(^ - - «S, 0} + (B.6) 

Po,j ' ^ 
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where the sign function is defined as follows: 

G M, sign(C) = 



1 if C > 
1 otherwise. 



Appendix B.3. Proximity operator of the temporal regularization function 

A simple expression of the proximity operator of function h is not avail- 
able. We thus propose to split this regularization term as a sum of two more 
tractable functions hi and /12: 

Nr/2 

t=l 

and 

Nr/2-1 

h2{0 = K J2 \\T*C^'-T*eX- (B-8) 

t=l 

Since hi (resp. /12) is separable w.r.t the time variable its proximity oper- 
ator can easily be calculated based on the proximity operator of each of the 
involved terms in the sum of Eq. f lB.7p (resp. Eq. ( IB.Sp ). 
Indeed, let us consider the following function 

^ : X ^ M (B.9) 



(B.IO) 



where = /^||T* • ||^ and H is the linear operator defined as 

H xC^ — ^ 

(a, b) a — b. 

Its associated adjoint operator i7* is therefore given by 

— ^C^ X (B.ll) 
a ^ (a, —a). 

Since we have HH* = 2Id, the proximity operator of ^ can easily be calcu- 
lated using (nl, Prop. 11]: 

prox^ = prox^^^ = Id + o (prox2^ - Id) o H. (B.12) 
The calculation of prox2^ is discussed in j9|. 
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Appendix B.^. Parallel Proximal Algorithm (PPXA) 
The function to be minimized has been reexpressed as 



:^st(C)=E 



E 



d\r) - S{r){TX'){r)\\l-. + g{C) 



t=l rG{l,...,X}x{l,...,y/i?}x{l,...,Z} 



+ /il(0 + /^2(C*)- 



(B.13) 



Since Jst is made up of more than two non-necessarily different iable terms, 
an appropriate solution for minimizing such an optimahty criterion is the 
PPXA [12]. In particular, it is important to note that this algorithm does 
not require subiterations as was the case for one of the algorithms proposed 
in j?!. In addition, the computations in this algorithm can be performed in a 
parallel manner and the convergence of the algorithm to an optimal solution 
to the minimization problem is guaranteed. 

The resulting algorithm for the minimization of the optimality criterion 
in Eq. (IB.lSp is given in Algorithm [H In this algorithm, the weights uji have 
been fixed to 1/4 for every i G {1, . . . , 4}. The parameter 7 has been set to 
200 since this value seems to give the fastest convergence in practice. The 
stopping parameter e has been set to 10~^. Using these parameters, the 
algorithm usually converges in less than 50 iterations. 
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